        subroutine nbodydec(rm)
c
c     Revision : 1.0
c
c input:  rm: Resonance mass
c
c     {\tt nbodydec} performs the decay of a resonance with mass rm in
c     its local rest  frame into nexit particles  with  4-momenta  and
c     masses  stored  in the array  pnew (see  comment to SUB jdecay).
c     The accessible many-body  phase-space is homogenously populated,
c     i.e. each configuration has equal probability. The theory behind
c     this  approach  can be  found in  M.M. Block and J.D. Jackson, Z.
c     Phys. C 3, 255 (1980). The original  routine is contained in  CPC
c     (Code ACGJ). It has been modified for uQMD purposes.
c     More documentation and better readability are to follow.

        implicit none
        integer j,i,imin1
        include 'newpart.f'
        real*8 rm
        real*8 p4loc(0:3), p1loc(3), ploc(3)
        real*8 M(mprt), MEFFloc(mprt), MASS
        real*8 M1,F1,M11


        real*8 z3, v2, ptot, costheta, z4, v1, reci1, z9, delx2, xi,
     +   z10, pp1dot,ener1,p2, p1s, z5, sintheta, phii,psin,
     +   p1sq, u, u1, ximin,ximax,delu,energy,pi,wmax,w,esys,z2,s,
     +   z8,reci,b, delxi,a,delz,ranf

        LOGICAL MASSLESS
        integer ntry

        p4loc(0) = rm
        p4loc(1) = 0d0
        p4loc(2) = 0d0
        p4loc(3) = 0d0
        wmax = 1d0

        PI=4d0*ATAN(1d0)
C
        M1=0d0 !Initialize M1=sum of all masses.
        MASSLESS=.FALSE.
C       Read masses.
        DO I=1,nexit
          m(i) = pnew(5,i)
          M1=M1+M(I) ! M1= sum of all masses.
        ENDDO
        IF (M1.EQ.0.) THEN  ! If massless.
          MASSLESS=.TRUE.
          WMAX=1d0
        ENDIF
        MEFFloc(1)=M(1) !Initialize Meff(1)
C       Initialize ESYS.
        ESYS=SQRT(p4loc(0)*p4loc(0)+p4loc(1)*p4loc(1)+
     +            p4loc(2)*p4loc(2)+p4loc(3)*p4loc(3))
C

C       Main Calculation
C
        ntry=0
60        W=1d0    !Initial weight, for each new event.
          ntry=ntry+1
          MASS=M1 !Initial M=total mass of ALL particles.
          ENERGY=p4loc(0)    !Initialize E to E*=cms energy.

          DO I=nexit-1,2,-1   !Loop over all N-2 effective masses needed.
            U1=M(I+1)/ENERGY
            U=U1**2
            MASS=MASS-M(I+1)
C           MASS=SUM of all rest masses of the REMAINING particles.
            XIMIN=(MASS/ENERGY)**2          !This is Xi,minimum.
            DELU=1d0-U1
            XIMAX=DELU**2           !This is Xi,maximum,where
                        !XIMIN <= XI <= XIMAX.
            DELXI=XIMAX-XIMIN       ! DELXI=delta (XI)=XIMAX-XIMIN.
            B=(1d0+U-XIMIN)          ! Used in FAST event generator.
            A=dble(I)*B            ! A=commonly used factor.
            IMIN1=I-1               ! IMIN1=commonly used factor.
            DELZ=(A-dble(IMIN1)*DELXI)*DELXI**IMIN1        ! DELZ=Zmax-Zmin
C
C           Here, we introduce the FAST generators.
C
            RECI=1d0/dble(I)
            S=1d0/(dble(I)-dble(IMIN1)*DELXI/B)    ! The probability
                        ! for the distribution i*(i-1)*y**(i-2)*(1-y).
100         Z2=ranf(0)
            IF (Z2.LT.S) THEN  ! The distribution i*(i-1)*y**(i-2)*(1-y).
                Z8=ranf(0)
                Z9=ranf(0)
                RECI1=1d0/dble(IMIN1)
                DELX2=DELXI*Z8**RECI*Z9**RECI1
              ELSE
                Z10=ranf(0)
                DELX2=DELXI*Z10**RECI   ! The probability distribution
                                ! i*y**(i-1)
            ENDIF
            IF (DELX2.EQ.0.) GOTO 100 ! Guards against division by 0.
ctp060202 110         XI=XIMIN+DELX2 ! XI=XImin+deltaXI.
            XI=XIMIN+DELX2      ! XI=XImin+deltaXI.
                        ! We reweight (multiply) W by
                                ! DELZ*F1*[(XI/(XI-XIMIN))**(I-2)]/(1+U-XI) .
            W=W*DELZ*F1(U,XI)*((XI/DELX2)**(I-2))/(1d0+U-XI)
            MEFFloc(I)=ENERGY*SQRT(XI) ! Store effective mass I, and
                                ! update E for next effective mass.
            ENERGY=MEFFloc(I)
          ENDDO
          V1=(M(1)/ENERGY)**2     ! Set up final weight, with particles 1 and 2.
          V2=(M(2)/ENERGY)**2
          W=W*F1(V1,V2) ! We now have the FINAL weight.
                !We find WMAX, the max weight, here.
c          IF (W.GT.WMAX) WMAX=W   ! Update WMAX.
                        ! This routine selects W=1 (unweighted events).
          Z3=ranf(0)
          IF (W.LT.WMAX*Z3.and.ntry.le.1000) THEN
             GOTO 60
          ENDIF
        ! We have accepted event, so see if we Lorentz transform it.

        M11=p4loc(0)
        p1loc(1)=p4loc(1)
        p1loc(2)=p4loc(2)
        p1loc(3)=p4loc(3)

        !Iterate over all blob masses, MEFF(I), where MEFF(1)=M(1),MEFF(N)=E*.
        DO 2500 I=nexit,2,-1
            ENERGY=.5d0*(M11+(M(I)**2-MEFFloc(I-1)**2)/M11)
            PTOT=SQRT(ENERGY**2-M(I)**2)
                !Find RANDOM cos(theta*)=COSTHETA, random PHI*=PHI
                ! SINTHETA=SIN(THETHA*)
            Z4=ranf(0)
            COSTHETA=2d0*Z4-1d0 ! -PI <= THETA* <= PI
            SINTHETA=SQRT(1d0-COSTHETA**2)
            Z5=ranf(0)
            PHII=2d0*PI*Z5   ! 0 <= PHI* <= 2*PI, random PHII
            PSIN=PTOT*SINTHETA !Commonly used combination.
                ! Calculate momentum compon. of particle I, ploc(k), k=1 to 3.
            ploc(1)=PSIN*COS(PHII) ! x-component.
            ploc(2)=PSIN*SIN(PHII) !y-component.
            ploc(3)=PTOT*COSTHETA ! z-component.
            P1SQ=p1loc(1)**2+p1loc(2)**2+p1loc(3)**2
            P1S=SQRT(P1SQ)
            ENER1=SQRT(P1SQ+M11**2)
                ! Calculate Plab(i) =
                        !P*(i) + betagamma(i)*
                        ! [Energy + betagamma(j).ploc(j)/(gamma+1)],
                                !where . means DOT product, i,j=x,y,z.
            PP1DOT=ploc(1)*p1loc(1)+ploc(2)*p1loc(2)+ploc(3)*p1loc(3) ! DOT product.
            A=(ENERGY+PP1DOT/M11/(1d0+ENER1/M11))/M11
                ! Plab=P1 for particle I;store in matrix OUT(K,I,3), update
                                ! new M11 and new p1loc()=ploc()-p1loc().
            P2=0d0
            DO J=1,3
                 ploc(J)=ploc(J)+A*p1loc(J) !Store new ploc().
                 P2=P2+ploc(J)*ploc(J) !Get square of ploc() vector.
                 p1loc(J)=p1loc(J)-ploc(J) !Update p1loc()
            ENDDO
            ENERGY=SQRT(P2+M(I)**2) !Store new ENERGY.
            M11=MEFFloc(I-1) ! Update M11.
c            WRITE (5,2600) M(I),ENERGY,ploc(1),ploc(2),ploc(3)
            pnew(5,i) = m(i)
            pnew(4,i) = sqrt(m(i)**2+ploc(1)**2+ploc(2)**2+ploc(3)**2)
            pnew(1,i) = ploc(1)
            pnew(2,i) = ploc(2)
            pnew(3,i) = ploc(3)
2500    ENDDO
c2600    FORMAT(0P,F7.4,2X,G13.7,5X,G13.7,3X,G13.7,3X,G13.7)
        P2=0d0 ! Do LAST particle here.
        DO J=1,3
            ploc(J)=p1loc(J) ! Store last ploc().
            P2=P2+ploc(J)*ploc(J)
        ENDDO
        ENERGY=SQRT(P2+M(1)**2) ! Store last ENERGY.
c               WRITE (5,2600) M(1),ENERGY,ploc(1),ploc(2),ploc(3)
c               WRITE (5,*)
            pnew(5,1) = m(1)
            pnew(4,1) = sqrt(m(1)**2+ploc(1)**2+ploc(2)**2+ploc(3)**2)
            pnew(1,1) = ploc(1)
            pnew(2,1) = ploc(2)
            pnew(3,1) = ploc(3)

        RETURN
        END
!-------------------------------------------------------------------------
        FUNCTION F1(V1,V2)
        ! Function F1(V1,V2)=SQR(1+(V1-V2)**2-2*(V1+V2))=2*(P*)/(E*).
        implicit none
        REAL*8 F1, F2, V1, V2
        F2=1d0+(V1-V2)**2-2d0*(V1+V2)
        IF (F2.LE.0d0) THEN
             F1=0d0
           ELSE
             F1=SQRT(F2) ! Guard against sqr(-).
        ENDIF
        END


      function M_inv_2(v01,vx1,vy1,vz1,
     +                 v02,vx2,vy2,vz2)
      real*8 M_inv_2,v01,vx1,vy1,vz1,
     +               v02,vx2,vy2,vz2

      M_inv_2 = sqrt((v01+v02)**2
     +              -(vx1+vx2)**2
     +              -(vy1+vy2)**2
     +              -(vz1+vz2)**2)
      return
      end

      function M_inv_3(v01,vx1,vy1,vz1,
     +                 v02,vx2,vy2,vz2,
     +                 v03,vx3,vy3,vz3)
      real*8 M_inv_3,v01,vx1,vy1,vz1,
     +               v02,vx2,vy2,vz2,
     +               v03,vx3,vy3,vz3

      M_inv_3 = sqrt((v01+v02+v03)**2
     +              -(vx1+vx2+vx3)**2
     +              -(vy1+vy2+vy3)**2
     +              -(vz1+vz2+vz3)**2)
      return
      end






      subroutine jdecay(rm)
C        input px,py,pz : CM-momenta of total system
C              rm:        Mass of resonance (sqrt(s))
c for pnew and pgen :
c      first index: 1=px, 2=py, 3=pz, 4=E, 5=m0
c      second index: particle number
      implicit none
       include 'newpart.f'
       real*8 pgen(5,mprt),rnd(mprt),u(3),beta(3),wt,tmp
       real*8 wtmax,rm,sum,pi,sum1,sum2,pcms,ranf,gamma,bp,phi,qcm,r1234

       parameter(pi=3.141592654)
       integer n,nadd1,i,j,ii,k
c
       pgen(1,1)=0d0
       pgen(2,1)=0d0
       pgen(3,1)=0d0
       pgen(5,1)=rm
       pgen(4,1)=rm
c
       nadd1=nexit-1
c
       pgen(5,nexit)=pnew(5,nexit)

c Two body decay
c ---------------
       if(nexit.eq.2) goto 400
       sum=0d0
c sum: sum of masses in the outgoing channel
       do 20 n=1,nexit
          sum=sum+pnew(5,n)
 20    continue

c     calculate maximum phase-space weight wtmax
c     ------------------------------------
       wtmax=0.5d0
       sum1=pgen(5,1)
       sum2=sum-pnew(5,1)
       do 200 i=1,nadd1
          wtmax=wtmax*pcms(sum1,sum2,pnew(5,i))
          sum1=sum1-pnew(5,i)
          sum2=sum2-pnew(5,i+1)
 200   continue

c     generate uniform nexit-body phase space
c     --------------------------------------
300   continue
c first generate nexit random numbers with decreasing value
c as excess energy distribution weights
      rnd(1)=ranf(1)
      do 110 i=2,nexit
         rnd(i)=ranf(1)
         do 120 j=i,2,-1
            if(rnd(j).gt.rnd(j-1)) then
               tmp=rnd(j-1)
               rnd(j-1)=rnd(j)
               rnd(j)=tmp
            endif
 120      continue
 110   continue
c last weight has to be zero
      rnd(nexit)=0d0
c now ?
      wt=1d0
      sum1=sum
      do 330 i=2,nexit
         sum1=sum1-pnew(5,i-1)
         pgen(5,i)=sum1+rnd(i)*(pgen(5,1)-sum)
         if(pgen(5,1)-sum.lt.0.0) write(6,*)'glrrrrrp'
         wt=wt*pcms(pgen(5,i-1),pgen(5,i),pnew(5,i-1))
 330  continue
      r1234=ranf(1)
      if(wt.lt.r1234*wtmax) goto 300

c     carry out two-body decays in pgen frames
c     ----------------------------------------
 400  continue
      do 410 i=1,nadd1
         qcm=pcms(pgen(5,i),pgen(5,i+1),pnew(5,i))
c        u(3) is cos(theta)
         u(3)=2d0*ranf(1)-1d0
         phi=2d0*pi*ranf(1)
         u(1)=sqrt(1d0-u(3)**2)*cos(phi)
         u(2)=sqrt(1d0-u(3)**2)*sin(phi)
         do 420 j=1,3
            pnew(j,i)=qcm*u(j)
            pgen(j,i+1)=-pnew(j,i)
 420     continue
         pnew(4,i)=sqrt(qcm**2+pnew(5,i)**2)
         pgen(4,i+1)=sqrt(qcm**2+pgen(5,i+1)**2)
 410  continue
      do 430 j=1,4
         pnew(j,nexit)=pgen(j,nexit)
 430  continue

c     boost pgen frames to lab frame
c     -------------------------------------------------
      do 500 ii=1,nadd1
         i=nexit-ii
         do 510 j=1,3
            beta(j)=pgen(j,i)/pgen(4,i)
 510     continue
         gamma=pgen(4,i)/pgen(5,i)
         do 520 k=i,nexit
            bp=beta(1)*pnew(1,k)+beta(2)*pnew(2,k)+beta(3)*pnew(3,k)
            do 530 j=1,3
               pnew(j,k)=pnew(j,k)+gamma*beta(j)*(pnew(4,k)
     &              +bp*gamma/(gamma+1d0))
 530        continue
            pnew(4,k)=gamma*(pnew(4,k)+bp)
 520     continue
 500  continue

      return
      end

